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The neutron superfluidity in the inner crust of a neutron star has been traditionally studied con- 
sidering either homogeneous neutron matter or only a small number of nucleons confined inside 
the spherical Wigner-Seitz cell. Drawing analogies with the recently discovered multi-band super- 
conductors, we have solved the anisotropic multi-band BCS gap equations with Bloch boundary 
conditions, thus providing a unified description taking consistently into account both the free neu- 
trons and the nuclear clusters. Calculations have been carried out using the effective interaction 
underlying our recent Hartree-Fock-Bogoliubov nuclear mass model HFB-16. We have found that 
even though the presence of inhomogeneities lowers the neutron pairing gaps, the reduction is much 
less than that predicted by previous calculations using the Wigner-Seitz approximation. We have 
studied the disappearance of superfluidity with increasing temperature. As an application we have 
calculated the neutron specific heat, which is an important ingredient for modeling the thermal 
evolution of newly-born neutron stars. This work provides a new scheme for realistic calculations 
of superfluidity in neutron-star crusts. 
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I. INTRODUCTION 

The possibility of superfluidity inside neutron stars was 
suggested a long time ago by Migdal [l|, only two years 
after the formulation of the theory of electron supercon- 
ductivity by Bardeen, Cooper and Schrieffer (BCS) [2J 
and before the discovery of the first pulsars. This predic- 
tion was later supported by the observation of the long 
relaxation time of the order of months following the first 
observed glitch in the Vela pulsar [3] . Glitches were sub- 
sequently observed in other pulsars. Pulsar glitches are 
believed to be related to the dynamics of the neutron 
superfluid permeating the inner layers of the solid neu- 
tron star crust Understanding the properties of 
this neutron superfluid is also of prime importance for 
modeling the cooling of newly-born neutron stars @-Qil 
and strongly magnetised neutron stars [Tlj j . the thermal 
relaxation of quasi-persistent X-ray transients [13, [l3| or 
the quasi-periodic oscillations recentl y d etected in the gi- 
ant flares of Soft-Gamma Repeaters |14l-tl6|. 

So far most microscopic studies of neutron superflu- 
idity have been devoted to the case of uniform infinite 
neutron matter [l7j . However the coherence length of 
the neutron superfluid in the neutron-star crust is typ- 
ically smaller than the lattice spacing and may even be 
comparable to the size of the nuclear clusters in some 
layers [H, [H (see also Section 8.2.3 of Ref. [20]). This 
situation is in sharp contrast to that encountered in or- 
dinary type I electron superconductors, for which the 
electron Cooper pairs are spatially extended over macro- 
scopic distances so that the order parameter is essentially 



uniform |2l| (note that neutron stars are much too hot 
for electrons to be superconducting there, see for instance 
the discussion in Section 8.1 of Ref. [13|). The effects of 
the inhomogeneities on neutron superfluidity in neutron- 
star crust have been studied in the mean-field approxi- 
mation with realistic nucleon-nucleon potential s 1221 and 
effective nucleon-nucleon interactions (1, [l(J [23|42d| . Sys- 
tematic fully self-consistent calculations in the entire in- 
ner crust (at zero temperature) have been recently car- 
ried out using semi- microscopic energy functionals [27| . 
In this latter work, it was found that even though inho- 
mogeneities are small in the densest regions of the crust, 
the 1 So neutron pairing gaps are strongly reduced com- 
pared to those in uniform neutron matter. Moreover the 
weaker the pairing force is, the stronger is the suppres- 
sion of the gaps. In somes cases, the gaps even almost 
completely vanish. All these quantum calculations have 
been carried out in the Wigner-Seitz (W-S) approxima- 
tion according to which the lattice is decomposed into a 
set of identical spherical cells centered around each clus- 
ter. The radius of each sphere is chosen so that its volume 
is equal to 1/pn, where p^ is the cluster density (num- 
ber of lattice sites per unit volume). As discussed in 
Ref. two different boundary conditions yielding an 
almost constant neutron density p n (r) near the cell edge, 
can be chosen. While the difference in the predicted pair- 
ing gaps are small for average nucleon densities p <j 0.03 
fm~ 3 , the uncertainties become increasingly large in the 
deeper layers of the crust [27j . This limitation of the W- 
S method is related to the existence of spurious neutron 
shell effects due to the discretisation of the single-particle 
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(s.p.) energy spectrum (28l.[29|. 

In this paper, we present the first calculations of neu- 
tron superfluidity in neutron-star crusts going beyond the 
W-S approach by using the BCS theory of anisotropic 
multi-band superconductivity. This theory is briefly re- 
viewed in Section|nl Our model of the neutron-star crust 
are discussed in Section Mil The solutions of the BCS 
equations are presented in Section IIVI We have focused 
on the deepest regions of the crust, for average nucleon 
densities p between 0.05 and 0.07 fm -3 , where the results 
from the W-S approximation are the most uncertain [27| . 
In Section [V] we discuss the validity of the local density 
approximation (LDA) and the importance of proximity 
effects. The disappearance of superfluidity with increas- 
ing temperature is studied in Sections IVII and IVIII In 
Section [Villi we present numerical results of the neutron 
specific heat. 

II. BCS THEORY OF SUPERFLUID 
NEUTRONS IN NEUTRON-STAR CRUSTS 

The standard formulation of the BCS theory starts 
with the Hamiltonian [2] 

CT,Q,fc 

+ X] V ak0k' c ifct C i-*4 C /3-fc'4- C /3fc't W 
a,f),k,kf 

where c^ fc (c akcr ) are the creation (annihilation) oper- 
ators for Bloch states with wave vector k, band index 
a and spin a, e ak are the s.p. energies (assumed to be 
independent of the spin state), p the chemical potential 
and V ak a k ' are the matrix elements of the two-body pair- 
ing interaction. In the mean-field approximation at finite 
temperature, the quasi-particle (q.p.) energies are given 
by 

E ak = ^{e ak -pf + ^ ak (2) 

where A ak are solutions of the anisotropic multi-band 
BCS gap equations (setting the Boltzmann constant 
fc B = 1) 

P k' P k 

In conventional superconductors the pairing interac- 
tion is induced by electron-phonon coupling 21]. It is 
usually a very good approximation to take the matrix 
elements V ak p k ' as constant and non-zero only within a 
small energy shell of the order ~ huj p around the Fermi 
level, where w p is the ion-plasma frequency. In this case 
the gap parameters A ak depend neither on the band in- 
dex a nor on the wave vector k and are all equal to a 
single constant, the pairing gap A Q. The possibility 



of multi-band superconductors characterized by the ex- 
istence of several pairing gaps A a was raised soon after 
the formulation of the BCS theory [3(|, but clear exper- 
imental evidence was lacking until the discovery in 2001 
of superconductivity in magnesium diboride, whose un- 
usual properties can be nicely explained by a two-band 
model [3l[ l32j. Since then many other multi-band su- 
perconductors have been found such as the iron pnictide 
superconductors (33|. In these materials, several bands 
can intersect the Fermi level, yielding a complex multi- 
sheeted Fermi surface. Pairing is still thought to be me- 
diated by the exchange of phonons but the electrons on 
the different sheets of the Fermi surface feel very different 
pairing interactions leading to the existence of different 
gaps. In neutron-star crust, the formation of neutron 
pairs giving rise to superfluidity is directly triggered by 
the strong neutron-neutron interaction which is always 
attractive at low densities in the 1 Sq channel. The num- 
ber of bands contributing appreciably to the pairing gap 
A ak can thus be huge (about ~ 10 2 — 10 3 in the dense 
layers of the inner crust considered in this work) . Since 
the matrix elements of the pairing force may a priori vary 
appreciably, we have solved Eqs. ([3]) in the most general 
case. 



III. MODEL OF NEUTRON-STAR CRUST 

We have determined the equilibrium structure and 
composition of the inner crust of neutron stars by using 
the fourth-order Extended Thomas-Fermi method with 
quantum shell effects added via the Strutinsky-Intcgral 
theorem. This so-called ETFSI method, as applied to 
the equation of state of neutron-star crusts, has been 
described in detail in Ref. 34]. It is a high-speed ap- 
proximation to the self-consistent Hartree-Fock method. 
We have neglected the small neutron shell effects [35j, 
the estimation of which in current calculations is plagued 
by the approximate treatment of the interaction between 
the unbound neutrons and the nuclear lattice [HI, [29| . 
The calculations have been carried out using an effective 
nucleon-nucleon interaction of the Skyrme type 

v S ^(ri,rj) = t (l + x P a )5{r tJ ) 

+ -t x {\ + x 1 P rJ )X 1 [p% 8{r t] ) + 5( rij )pl] 

+ t 2 (l + x 2 P<T)-oPij.S(rij)pij 
n 

+ h 3 (l + x 3 P a )p(ryS( rij ) 

T2W Q (o-i + (Tj) -Pij X 5(rij)pij , (4) 

n 

where = n-rj, r = (ri+rj)/2,Pij = -i^(Vi-Vj)/2 
is the relative momentum, and P a is the two-body spin- 
exchange operator and p(r) is the total nucleon density 
at position r. The pairing interaction that we take here 
acts only between nucleons of the same charge state q 
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(q = n or p for neutron or proton, respectively) and is 
given by 

v^(n, rj ) = v*«[p n (r),p p (r)} S( rij ) , (5) 

where m = ri—rj, r = (r»+rj)/2 and p n (r) and p p (r) are 
the neutron and proton density at position r respectively. 

We have adopted the parametrisation BSkl6, underly- 
ing the HFB-16 nuclear mass model 36]. The parame- 
ters of this force are given in Table |U This force is par- 
ticularly suitable for studying neutron-rich environments 
such as neutron-star crusts since it has been constrained 
to reproduce the equation of state and the 1 S'o pairing 
gap of infinite homogeneous neutron matter, as calcu- 
lated for the realistic Argonne V14 potential and shown in 
Fig. [TJ The expression of v nq [p n (r), p p (r)] can be found 



in Ref. [36[. Moreover, we can hope that the nuclear in- 
homogeneities in the neutron-star crust will be properly 
taken into account, given the excellent fit to essentially 
all the available experimental nuclear mass data. 

In order to solve the BCS Eqs. ([3]), we first need to de- 
termine the neutron s.p. energies e ak and wavefunctions 
i Pak( r )- F° r this purpose, we have solved the following 
three-dimensional Schroedinger equation 

h 2 

-V- 2M «/ f x Vy afc (r) + Z7 n (r)y afc (r) = e ak (p ak (r) . (6) 

The effective mass M*(r) and the potential U n (r) are 
given by (we neglect the small rearrangement term com- 
ing from the pairing force, see for instance 
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where r n (r) and r p (r) are the neutron and proton kinetic- 
energy densities respectively (we have introduced r = 
T n + t p ). Both M*(r) and U n (r) were determined using 
the ETFSI fields. In Eq. ©, we have neglected the spin- 
orbit coupling arising from the last term of Eq. ((4]) . This 
approximation is justified because the spin-orbit coupling 
is proportional to Vp„ and Vp P (see for instance Ap- 
pendix A of Ref. [36(). In neutron-star crust, nuclear 
clusters have a very diffuse surface as can be seen in Fig|2] 
and consequently the spin-orbit coupling is much smaller 
than that in isolated nuclei {29|. 

Although neutron shell effects represent a small correc- 
tion to the total energy density [35|, they are expected 
to have a much stronger impact on neutron superfluidity 
due to the highly non-linear nature of the pairing phe- 
nomenon. This is the reason why we have not followed 
the usual practice of applying the W-S approximation 
for solving Eq. (|6]), but have imposed the Bloch bound- 
ary conditions 

<Pak(r + e) = exp(i& ■£)>p a k(r) , (9) 

where £ denotes any lattice translation vector. Note that 
in this case, Eq. © has to be solved for each wave vector 



k while in the W-S method only a single wave vector is 
considered, namely k = (see Ref. [29| for a discussion 
about the W-S approximation). Following the standard 
assumptions, we have considered a body-centered cubic 
lattice The BCS Hamiltonian ((TJ can then be ob- 

tained from the s.p. states, once the pairing interaction 
has been specified. The matrix elements of the pairing 
force (JSJ between Bloch states are given by an integral 
over the W-S cell of volume V ce ii 

V ak0k '= f d 3 rv*[p n (r)}\v ak (r)\ 2 \^ k >(r)\ 2 , (10) 
Jws 

the Bloch wavefunctions ip ak (r) being normalized accord- 
ing to 

/ d 3 r|^ Qfc (r)| 2 = l. (11) 
Jws 

The W-S cell that we consider here is a truncated octa- 
hedron, as determined by the body-centered cubic lattice 
geometry. It should not be confused with the spherical 
cell used in the W-S approximation (29j . 
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TABLE I: Skyrme parameters of the force BSkl6 [36j. £a is 
a s.p. energy cutoff above the chemical potential introduced 
in order to regularize the divergences associated with the zero 
range of the pairing force. 



t [McV fm 3 ] 


-1837.23 


ti [McV fm 5 ] 


383.521 


tl [McV fm 5 ] 


-3.41736 


t 3 [McV fm 3+3 T] 


11523.0 


X 


0.432600 


Xl 


-0.824106 


X2 


44.6520 


X3 


0.689797 


Wo [McV fm 5 ] 


141.100 


7 


0.3 


£A [MeV] 


16 




kp [fm -1 ] 

FIG. 1: So neutron pairing gap in infinite uniform neutron 
matter vs Fermi wave number k-p = (37r 2 p n ) 1/ ^ 3 , as used in 
model HFB-16 [H. 



IV. NEUTRON PAIRING GAPS 

We have considered five different layers of the inner 
crust in the average nucleon density range between 0.05 
fm~ 3 and 0.07 fm" 3 . Results of the ETFSI calculations 
at T = arc summarized in Table |TT] and the nucleon 
density profiles are plotted in Fig. [5J 

We have solved iteratively the anisotropic multi-band 
BCS gap equations (J3J) for each average nucleon density 
p. Results are summarized in Table HTT1 Due to the ultra- 
violet divergence induced by the zero range of the pair- 
ing interaction ([5]), the summation in Eq. Q has to be 
truncated. We have imposed the same s.p. energy cutoff 
e\ = 16 MeV above the chemical potential as used in the 
determination of the BSkl6 force through optimization of 
the mass fit. For each temperature the chemical potential 
has been recalculated neglecting pairing since the pairing 
gaps are much smaller than the Fermi e nerg y (note that 
the same approximation was used in Ref. |36| to construct 



TABLE II: Ground-state composition of the neutron-star 
crust using the ETFSI method with Skyrme force BSkl6. 
p is the average nucleon density, Z and A the equilibrium 
numbers of protons and nucleons in the W-S cell respectively, 
and f>Bn the neutron background density outside clusters (sec 
Ref. HI). 



P [fm- 3 ] 


Z 


A 


p Bn [mi -3 ] 


0.070 


40 


1258 


0.060 


0.065 


40 


1264 


0.056 


0.060 


40 


1260 


0.051 


0.055 


40 


1294 


0.047 


0.050 


40 


1304 


0.043 



the effective density-dependent pairing strength from the 
pairing gap in infinite homogeneous neutron matter) . We 
have solved Eq. © by expanding the s.p. wavefunctions 
into plane-waves 

<Pak{r) =cxp(ifc-r)^(^ Qfe (G)exp(iG-r) (12) 

G 

in which G are reciprocal lattice vectors. Since by defi- 
nition 

exp(iG-f) = l (13) 

for any vectors G and £, the Bloch boundary condi- 
tions ([SJl are automatically satisfied. We have included 
all Fourier components with reciprocal lattice vectors G 
such that \k + G\ < Q. Q has been adjusted so that the 
s.p. energies are computed with an accuracy of a few 
keV. We have evaluated the summation in Eq. Q us- 
ing the special-point method (3?J ■ We have also applied 
this method to compute the pairing matrix elements (jlOl) . 
On general grounds one may expect that \f a k( r )\ 2 , and 
thereby V ak/3k ' and A akl are weakly dependent on k, 
since bound states are vanishingly small outside clus- 
ters where Bloch boundary conditions are imposed while 
continuum states depend on k essentially through only a 
phase factor exp(iA; ■ r). Indeed we have found that the 
summation in Eq. (J3]) converges quickly with the num- 
ber of special fc-points. An error below 1% for the av- 
eraged pairing gap Ap can be reached with 30 Appoints. 
Keeping just one term, corresponding to the mean- value 
point [38[ , yields a result with a few percent precision, as 
can be seen in Table IIIII This method has indeed proved 
to be surprisingly accurate in solid state physics to com- 
pute the electron density and dielectric matrix [39( . The 
convergence of the real-space integrations in Eq. (1101) is 
slower because of the oscillating behavior of the wave- 
functions. We have checked that the solutions of the 
BCS gap equations converge to a few keV accuracy using 
110 special r-points. 

As can be seen in Fig. [31 the dependence of the pair- 
ing gaps A a k on the band index a and wave vector k 
is quite significant. At the Fermi level the pairing gaps 
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FIG. 2: Equilibrium neutron number density p n (r) (solid line) and proton number density p P (r) (dashed line) inside the 
Wigner-Seitz cell of different neutron-star crust layers with average nucleon density J>. The densities, expressed in fm -3 , have 
been obtained from the ETFSI method Q at T = with the Skyrme force BSkl6 Q^. 



vary by about ~ 0.2 — 0.4 MeV. While all s.p. states lying 
in the continuum contribute to the average gap Ap, we 
have found that Ap remains almost unchanged if bound 
states are excluded from the summation in Eq. ((3|) . How- 
ever this does not imply that neutrons inside clusters do 
not have any impact on the average pairing gap. The 
inhomogeneous distribution of neutrons modifies the s.p. 
energies and the matrix elements of the pairing force. 



For comparison, we have calculated the neutron pair- 
ing gap A u without nuclear clusters, assuming that un- 
bound neutrons are uniformly distributed with the den- 
sity pBn corresponding to the neutron background den- 
sity outside clusters (values of pBn for the different 
crustal layers are indicated in Table [TTJ) . In this limit- 
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2 
1.8 



> 1 6 



-a 1.4 



(a) 



1 1 ^^d^WWRwRJ 




40 -30 -20 -10 10 20 




-40 -30 -20 -10 10 20 
e ak [MeV] 



40 -30 -20 -10 10 20 
e ak [MeV] 




e ak [MeV] 

FIG. 3: Neutron pairing gaps A a k vs s.p. energies e a k for the different crustal layers. The arrow indicates the position of the 
chemical potential at T = 0. The BCS equations ([3| have been solved at T — together with Eq. using the neutron-star 
crust composition shown in Fig. [2] The dashed line and the dotted line represent the gaps A u and A u respectively given in 
Table Hill 



ing case, the Bloch wavefunctions reduce to plane waves 



<Pak (r) 



Vo 



■exp[i(k + G a ) -r] , (14) 



trix elements become all equal 

v*[p B n] 



V, 



ak0k' 



V C ell 



(15) 



where G a are reciprocal lattice vectors. The pairing ma- 



As a result, the pairing gaps are independent of a and 
k, and are the solutions of the usual isotropic BCS equa- 



7 



TABLE III: Neutron pairing gaps in the neutron star crust 
for different average nucleon density p at T — 0. Ap (Apo) 
is the pairing gap obtained from Eq. (|3]) after averaging over 
continuum states and using the special point method (mean 
value point method) for the summation over k. A u and A u 
are the pairing gaps in uniform infinite neutron matter for 
the neutron density pg„ and for the average neutron density 
Pn = (A— Z)/VoeU respectively. 



p [fm- s ] 


A F [MeV] 


Apo [MeV] 


A u [MeV] 


A u [MeV] 


0.070 


1.44 


1.39 


1.79 


1.43 


0.065 


1.65 


1.59 


1.99 


1.65 


0.060 


1.86 


1.81 


2.20 


1.87 


0.055 


2.08 


2.07 


2.40 


2.10 


0.050 


2.29 


2.27 


2.59 


2.33 



tions 



V*\PL 



8tt 2 



2M 



Bn 



3/2 



, d£ W) 



tanh 



with 



E(e) = V(e-e BF ) 2 



A 2 



in which 



9 M* 



E(e) 
2T ' 
(16) 



(17) 



(18) 



is the Fermi energy, fcps = (37r 2 psn) 1 ^ 3 and M Bn is the 
effective mass in neutron matter at density pBn- From 
the definition of the pairing strength, it follows immedi- 
ately that A u at T = 0, is nothing else but the micro- 
scopic neutron pairing gap, shown in Fig. [I] evaluated at 
the neutron density p n = psn- Results for the different 
crust layers are shown in Table IIIII Comparing Ap and 
A u , it can be seen that the presence of inhomogeneities 
lowers the averaged neutron pairing gap by ~ 10 — 20%. 
This reduction is much smaller than that found in previ- 
ous calculations based on the W-S approach [27]. For 
instance, for p = 0.058 fm -3 , Ap was suppressed by 
~ 40 — 60% depending on the choice of boundary condi- 
tions (last line of Table 2 in Ref. [13]). 



V. PAIRING FIELD OF THE NEUTRON 
SUPERFLUID AT ZERO TEMPERATURE AND 
LOCAL DENSITY APPROXIMATION 

The effects of the inhomogeneities on the neutron su- 
perfluid can be more directly seen by computing the neu- 
tron pairing field, defined by 



A n (r) 



f 



v™\p n {r),pp{r)}p n {r), 



(19) 



where p n (r) and p n (r) are the local normal neutron den- 
sity and abnormal neutron density respectively given (at 
T = 0)by 



1 



£ak — M 



Pn(r) 



A 

E 



<Pak(r)\ -—, 

E a k 



(20) 



(21) 



where the superscript A is to indicate that the summation 
includes only states whose s.p. energy lies below p + s\. 
We have computed the angle-averaged pairing field inside 
each W-S cell, as given by 

A„(r)=/ ^A„(r). (22) 

For each value of the radial coordinate r, we have per- 
formed the solid-angle integration using 30 uniformly dis- 
tributed spiral points on the sphere of radius r [4(J . In 
order to minimize the amount of computations, we have 
applied the mean- value point method for solving the BCS 
equations and for calculating the normal and abnormal 
densities. For comparison, we have also calculated the 
pairing field in the LDA, i.e., assuming that at each point 
r the pairing field is locally the same as that in uniform 
neutron matter for the density p n { r ) 

(23) 

As can be seen in Fig. 21 the LDA overestimates the spa- 
tial dependence of the pairing field thus indicating that 
pairing is highly non-local. Indeed the local coherence 
length of the neutron superftuid, defined by 



Al LDA )(r) = A u (p„(r)). 



i{r) = 



?i 2 k F (r) 



7rM*(r)Ai LDA) (r) 



(24) 



where 

fc F (r) = (3^ 2 pn(r)) 1/3 , (25) 

is larger than the size of the clusters, as shown in Fig. [5] 
This means that even if the center of mass of a Cooper 
pair is located outside clusters, one of the partner may ac- 
tually lie inside so that all neutrons are actually involved 
in the pairing process. As a result, pairing correlations 
are strongly enhanced inside clusters but are reduced in 
the intersticial region, leading to a smooth spatial vari- 
ation of the pairing field. These so-called proximity ef- 
fects are the most spectacular in the shallowest layer at 
p = 0.05 fm~ 3 where the neutron pairing field is increased 
by a factor ~ 4 inside clusters. The entire inner crust is 
therefore permeated by the neutron superffuid, including 
in the region occupied by clusters themselves. 

The previous considerations suggest to estimate Ap by 
calculating the pairing gap A u in uniform neutron matter 
for the average neutron density p n = N/V C ett instead of 
the neutron background density pBn- As can be seen in 
Table IIIII this gap indeed provides a very good approxi- 
mation to the average gap calculated numerically. 
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FIG. 4: Angle-averaged neutron pairing field A(r) (solid line) and LDA pairing field (dashed line) for the different crustal 
layers at T = 0. 



VI. CRITICAL TEMPERATURE 



represented by [4l| 



While the pairing gaps obtained from Eqs. ([3]) vary 
from one band to the other as shown in Fig. [31 we have 
found that they all share the same universal temper- 
ature dependence which (for pairing) can be well- 



A afc (T<T c )~ A afc (0)Wl- 



(26) 



with S ~ 3.23 and the critical temperature T c , is de- 
fined by the condition A a k(T > T c ) = for all states. 
For simplicity, we have assumed that the composition re- 
mains unchanged at finite temperature so that the s.p. 
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energies e a k and the matrix elements V ak p k > are inde- 
pendent of T. We have determined T c for each crustal 
layer by solving numerically Eqs. ([3]). For comparison, 
we have also calculated the critical temperature T cu of a 
uniform neutron superfluid with the corresponding neu- 
tron background density pg n . As indicated in Table HVl 
the actual critical temperature T c is systematically lower 
than T cu , as could have been expected from the results 
on the pairing gaps discussed in Section [Vj The critical 



temperature T cu of a uniform neutron superfluid with the 
average neutron density p n is much closer to T c . 

It is well-known from the BCS theory of superconduc- 
tivity that the ratio of the critical temperature to the 
pairing gap is universal [2J. We have checked that this 
relation holds in uniform neutron matter, namely 

T CU = ^A U , (27) 

(likewise for T cu ) where ( ~ 0.577 is the Euler- 
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Mascheroni constant (this result is not specific to the 
pairing model we have used. It has also been found in 
microscopic calculations using realistic nucleon-nucleon 
interactions, as discussed in Ref. 42]). What is perhaps 
more surprising is that we have found the same relation 
for the inhomogeneous neutron superfluid in the neutron- 
star crust 



exp(C) 



(28) 



(note that in the original BCS theory Q, the supercon- 
ductor is supposed to be isotropic and uniform). Since 
the band theory includes the limiting case in which all 
neutrons are bound inside clusters, we expect Eq. ([28]) to 
remain valid in finite nuclei. Quite interestingly, this con- 
clusion seems to be supported by self-consistent mean- 
field calculations in tin isotopes |43| even though in this 
case the averaged gap Af has to be suitably defined. 

We can understand Eq. (|2"5)) by considering a simple 
pairing model. Since Eq. (|28p holds in very different 
situations, this means that unlike the state-dependent 
pairing gaps A ak , the critical temperature is rather in- 
sensitive to the precise nature of s.p. states (see also 
Ref. [44| for a mathematical discussion in the context 
of multi-band superconductivity). This suggests replac- 
ing the pairing matrix V ak p k > by some constant coupling 
(V)/V ce ll- In this case the pairing gaps A ak become all 
equal to the same value A and Eq. ([3]) reduces to the 
finite-temperature isotropic BCS gap equations 



1 = _i.W 



2 V c 



de- 



al Jo 



9(e) 
E(e) 



tanh 



E(e) 
2T 



(29) 



where E(e) = y/ (e — fi) 2 + A 2 is the q.p. energy and 
g(e) is the density of s.p. states (for a given spin state) 
defined by 



9(e) =^<5(e Qfc -e)- 



(30) 



TABLE IV: Critical temperature for the onset of neutron su- 
perfluidity in the neutron-star crust for different average nu- 
cleon density p using different approximations: T c is the crit- 
ical temperature obtained after solving numerically the BCS 
Eqs. ([3]) while T cu and T cu are the critical temperatures of a 
uniform superfluid with density pBn and p„ respectively. All 
temperatures are indicated in units 10 10 K. 



p [far 3 ] 


T 

cu 


T c 


T 

J- cu 


0.070 


1.17 


0.96 


0.94 


0.065 


1.31 


1.09 


1.08 


0.060 


1.45 


1.24 


1.23 


0.055 


1.58 


1.38 


1.38 


0.050 


1.70 


1.52 


1.53 



with k F = (37r 2 p„) 1/3 . Since T c < e F , we take /i(T c ) ~ 
£f remembering that the lowest order correction is only 
of order (T c /e F ) 2 . Since the most important contribution 
to the integral in Eq. (I2U1) comes from s.p. states lying 
in the vicinity of the Fermi surface, we replace g(e) by 
g{e F ). For T = T c , A = and Eq. $ZS§ thus becomes 



An 2 h 



2*2 



(V}M*h 



£a/2T c 



d.r 



tanh I 



ef/2T 



(33) 



After remarking that T c <C £a, the integral in Eq. ([33 
can be solved analytically using 



y , tanhx , /4w\ 
dx ~ log ( -f ) + C 

o x 



7T 



(34) 



for y > 1. Replacing Eq. ([34]) in Eq. (|33|l . we find that 
the critical temperature is given by 



2exp(0 . 

T c = V £ F£Aexp 



87T 2 ?i 



2*2 



k F M*(V) 



(35) 



The critical temperature T c is mainly determined by 
unbound neutron s.p. states (see the discussion in Sec- 
tion [Vj. However, as shown in a previous work [4o| . the 
density of unbound neutron s.p. states is essentially un- 
affected by the inhomogencitics on an energy scale larger 
than a few hundred keV. Since the pairing gaps are typ- 
ically of the order of MeV, we can replace g(e) by its 
expression in uniform neutron matter with density p n 



9(e) 



Veen /2M* 



4tt 2 \ h 2 



3/2 



(31) 



where M* is the corresponding neutron effective mass. 
At T — the neutron chemical potential /i„ is given by 
the neutron Fermi energy 



£ F 



n 2 4 

2M* 



(32) 



On the other hand, we assume that the average pairing 
gap Ap can be obtained from the solution of the BCS gap 
Eq. ((3]) at T = after substituting the pairing matrix 
V ak p k i by the same pairing constant (V)/V ce n that was 
introduced previously for calculating T c . Eq. (3) thus 
reduces to Eq. ([29| which, at T = 0, reads 



2fi2 



(V)M*k F 



ea/A 



da^l+x 2 )- 1 / 2 , (36) 



where we have taken the density of states out of the in- 
tegral. Solving Eq. ([36]) for the pairing gap yields 



Ap = A = 2^/e F e\ exp 



2*2 



in 2 h 



k F M*(V) 

Comparing Eqs. (|3"5"|) and (|37p leads to Eq. 



(37) 
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VII. FINITE TEMPERATURE EFFECTS ON 
THE NEUTRON PAIRING FIELD 

At finite temperatures, the pairing field is still given by 
the same expression (ITS)) as for T = 0. But the normal 
neutron density and abnormal neutron density are now 
given by 



a.k 



E n 



Aanh% 
IT 



Pn(r) = V |^fc(r)| 2 ^- tanh E ' 

, E a k 

a.k 



ak 

2T 



(38) 



(39) 



In the LDA, the corresponding pairing field is given by 



A( LDA >(r,T) = A u (p„(r)Wl 



T 



T cu (r) 



for T < T cu (r) with 



exp(C) 



while for T > T cu {r), Ai LDA) (r,T) = 0. Fig. © shows 
the pairing field for different temperatures. The tem- 
perature dependence of the pairing field comes mainly 
from that of the abnormal density. As the tempera- 
ture gets closer to the critical temperature T c , the su- 
perfluid becomes more and more homogeneous. This is 
because the coherence length obtained after substitut- 
ing Al LDA) (r) by Ai LDA) (r,T) in Eq (JUJ), increases with 
T and even diverges when the local critical temperature 
is reached. This means that the LDA becomes worse 
at finite temperatures, strongly overestimating the im- 
pact of inhomogeneities. In particular Eq (1401) implies 

that Ai LDA ^(r,T cu (0)) vanishes inside clusters while it 
remains finite outside, leading to a sharp variation of the 
pairing field at the cluster surface. 



A u(Pn(r)) : 



(40) 



(41) 



VIII. SPECIFIC HEAT OF SUPERFLUID 
NEUTRONS 

Observations of the thermal X-ray emission from 
newly-born isolated neutron stars can potentially pro- 
vide valuable information on the structure of neutron- 
star crusts. Due to its relatively low neutrino emissivity, 
the crust of the neutron star cools less rapidly than the 
core and thus stays hotter. As a result, the surface tem- 
perature decreases slowly during the first ten to hundred 
years after the formation of the neutron star in a super- 
nova explosion and then drops sharply when the cooling 
wave from the core reaches the surface @, Q • The oc- 
currence of neutron superfluidity in the inner crust has a 



strong influence on the neutron specific heat and in turn, 
on the thermal relaxation time of the crust [t| [l(| • 



In Ref. [45j, we have shown that the specific heat of 
normal neutrons in the outermost layers of the inner crust 
is almost unaffected by the nuclear lattice and is essen- 
tially given by the specific heat of a uniform gas. At 
low temperatures T < T c , pairing correlations have to be 
taken into account. This is usually done by introducing 
a multiplication factor 1Z, defined by 



where Cy and Cy' are the specific heat of superfluid 
and normal neutrons respectively. Assuming that the free 
neutrons are uniformly distributed, the factor 7t is well- 
approximated by the analytical expression of Ref. [46[ 
(used in neutron-star cooling simulations) 



C ( y Sn) = KC^ 



(42) 



0.4186 + Vl-007 2 + (0.501u) 2 
x exp(1.456 - \JlAb& 2 + u 2 ) (43) 



with y = T/T cu and 



yjl-y[ 1.456 



0.157 1.764 



Vv 



(44) 



In order to assess the validity of these expressions for the 
inhomogeneous matter of neutron-star crusts, we have 
evaluated the neutron specific heat from the numerical 
solutions of the anisotropic multi-band BCS gap equa- 
tions ([3]) . We have applied the following expression from 
Ref. 47] using the temperature dependence of the gaps 
given by Eq. (|26|) 



1 \ - exp(E ak /T) 
Vcdi [l + exp^/T)] 2 

_f__d_ 

~ 2TdT 



Egk 

T 



A Qfe (T) : 



(45) 



Numerical results are shown in Fig. [7J The summation 
in ft-space in Eq. (|45l) has been carried out using the 
special-point method [37J. 

At low temperatures T -C T c , neutron pairing correla- 
tions are very strong. The pairing gaps are weakly de- 
pendent on T and the last term of Eq. (|45j) is negligible 
(de Genne approximation [47j)- In this regime, the neu- 
tron specific heat is exponentially suppressed compared 
to that of normal neutrons and decays approximately like 
exp(— Ap/T). However as the temperature gets higher, 
the contribution of the last term of Eq. ([4"5"j) becomes 
increasingly large, eventually leading to a sharp rise of 
the specific heat for T < T c . For T > T c , all pairing 
gaps vanish and neutrons superfluidity is destroyed. For 
T > T c and T <C £f, the specific heat increases linearly 
with T as can be clearly seen in Fig. [7] As in our previous 
work W5ji we have found that the numerical results for 
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FIG. 6: Angle-averaged neutron pairing field A(r) for the different crustal layers at different temperatures (T is indicated above 
each curve in 10 10 K). 



the normal neutron specific heat are in good agreement 
with the expression for a uniform gas 



r<( n ) — 



3h 2 



-T. 



(46) 



Note however that for the dense crustal layers consid- 
ered here, using M* and kp instead of Mg n and fcpB 
(as in Ref. (45|) yields a slightly better fit. As shown 
in Fig. [71 the transition between the superfluid and the 



normal regime is very sharp. The specific heat exhibits 
a discontinuity at T = T Cl approximately given by 



4 ra) (T c )-C 



( v\T c ) 



4")(T C ) 



\& exp(- 



2C) 



(47) 



after replacing A a f. by Ap in Eq. (T43t and the density of 
s.p. state by Eq. ([51]). Note that for T < T c , the neutron 
specific heat is enhanced by pairing, i.e. 1Z > 1. 
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FIG. 7: Specific heat of superfluid neutrons in different layers of the inner crust of a neutron star. The solid line has been 
obtained from Eq. (|45|) after solving numerically the BCS Eqs. ((3J The dashed and dotted lines correspond to homogeneous 
neutron matter with two different values of the critical temperature T c and T cu , respectively. 



We have found that Cy is numerically close to the 
expression for a uniform superfluid, namely Eqs. (142ft . 
P3"]) and (|4"fjj) . provided that the critical temperature is 
suitably renormalized (i.e., T cu must be replaced by T c ) 
as shown in Fig. [7] (renormalized homogeneous superfluid 
curve). 



IX. CONCLUSIONS 

We have clarified the effects of the nuclear inhomo- 
geneities on the neutron superfluidity in the deep layers 
of the neutron-star crust by solving the anisotropic multi- 
band BCS gap equations ([3]). We have properly taken 
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into account the interactions between the neutron super- 
fluid and the nuclear crystal lattice by imposing Bloch 
boundary conditions instead of using the W-S approxi- 
mation. Due to the presence of the nuclear clusters, neu- 
trons belonging to different bands and having different 
Bloch wave vectors feel different pairing interactions thus 
leading to a dispersion of the neutron pairing gaps A Q £ of 
about ~ 0.2 — 0.4 MeV around the Fermi level, as shown 
in Fig. [3] The neutron pairing gap Ap averaged over all 
continuum states is reduced due to the presence of in- 
homogeneities, but much less than predicted by previous 
calculations based on the W-S approximation 27]. Un- 
like the individual gaps A Q £, Ap is essentially unaffected 
by the band structure and is very well approximated by 
the gap A u in uniform neutron matter using the aver- 
age neutron density p n . The reason for this agreement 
lies in the highly non-local character of the pairing phe- 
nomenon involving both bound and unbound neutrons, 
as revealed by numerical calculations of the neutron pair- 
ing field A„(r). As a consequence the LDA strongly over- 
estimates the spatial variations of the pairing field. The 
discrepancies are particularly large inside clusters where 
the LDA incorrectly predicts a quenching of pairing cor- 
relations. 

Solving the BCS equations at finite temperatures, we 
have found that the temperature dependence of the 
pairing gaps A a k is universal and well-represented by 
Eq. p6|. Unlike what one might have naively expected, 
the critical temperature T c is not determined by the 
largest state-dependent pairing gap around the Fermi 
level, but rather by the average pairing gap Ap. More- 
over the ratio T c /Ap is approximately given by the BCS 
value ~ 0.58 [2j. We have explained these results by solv- 
ing the isotropic BCS Eqs. (f29|) in the weak-coupling ap- 
proximation leading to the analytical expression (|55|) for 
the critical temperature. We have computed the neutron 



pairing field at finite temperatures and we have shown 
that the LDA becomes worse as the temperature is in- 
creased. 

We have computed the neutron specific heat, which is 
an important ingredient for modeling the thermal evolu- 
tion of newly-born neutron stars [7H10j. We have found 
that the specific heat is modified by band effects but it 
can be easily estimated from the expression in uniform 
neutron matter, namely Eqs. (|42j) . (|43[) and (|46[) . by sim- 
ply renormalising T c . 

The conclusions of the present work may change in 
shallower regions of the crust where the matter is more 
inhomogeneous, as suggested by a recent study in dirty 
superconductors [48j |. In addition we have neglected the 
change in composition with increasing temperature [34| 
which could affect the values of the pairing gaps (hence 
also the critical temperature) for T > 10 10 K. In par- 
ticular we expect to find deviations of the ratio T c /Ap 
from the BCS value~ 0.58 when these thermal effects are 
taken into account. We have also left aside the modifica- 
tions of the pairing gaps due to many-body effects beyond 
the BCS approximation [49j . These various issues will be 
addressed in future works using the multi-band approach 
presented in this paper. 
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